AGBRESA Cognition

Required Packages

library(knitr) # for kable
library(readxl) # read xls
library(DT) # Print table nicely in markdown
library(reshape2) # wide to long
library(plyr) # ldply
library(tidyverse) # for fct_relevel

library(lme4) # for adding mixed model functions to plot
library(lmerTest) # for adding mixed model functions to plot
library(emmeans) # for adding mixed model marignal means and se to plot
library(qdapTools) # list_df2df(tbl_mm_fac_aov)

library(ggpubr) # for compare means, descriptives and theme_pubr
library(ggplot2)
library(ggsci)
library(scales)
library(cowplot) # for plot_grid, combine plots using cowplot

Print Session Info

print(sessionInfo())

## R version 4.1.1 (2021-08-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] cowplot_1.1.1   scales_1.1.1    ggsci_2.9       ggpubr_0.4.0   
##  [5] qdapTools_1.3.5 emmeans_1.7.0   lmerTest_3.1-3  lme4_1.1-27.1  
##  [9] Matrix_1.3-4    forcats_0.5.1   stringr_1.4.0   dplyr_1.0.7    
## [13] purrr_0.3.4     readr_2.0.2     tidyr_1.1.4     tibble_3.1.5   
## [17] ggplot2_3.3.5   tidyverse_1.3.1 plyr_1.8.6      reshape2_1.4.4 
## [21] DT_0.19         readxl_1.3.1    knitr_1.36     
## 
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-152        bitops_1.0-7        fs_1.5.0           
##  [4] lubridate_1.8.0     httr_1.4.2          numDeriv_2016.8-1.1
##  [7] tools_4.1.1         backports_1.2.1     utf8_1.2.2         
## [10] R6_2.5.1            DBI_1.1.1           colorspace_2.0-2   
## [13] withr_2.4.2         tidyselect_1.1.1    curl_4.3.2         
## [16] compiler_4.1.1      chron_2.3-56        cli_3.0.1          
## [19] rvest_1.0.2         xml2_1.3.2          mvtnorm_1.1-3      
## [22] digest_0.6.28       foreign_0.8-81      minqa_1.2.4        
## [25] rmarkdown_2.11      rio_0.5.27          pkgconfig_2.0.3    
## [28] htmltools_0.5.2     dbplyr_2.1.1        fastmap_1.1.0      
## [31] htmlwidgets_1.5.4   rlang_0.4.12        rstudioapi_0.13    
## [34] jquerylib_0.1.4     generics_0.1.0      jsonlite_1.7.2     
## [37] zip_2.2.0           car_3.0-11          RCurl_1.98-1.5     
## [40] magrittr_2.0.1      Rcpp_1.0.7          munsell_0.5.0      
## [43] fansi_0.5.0         abind_1.4-5         lifecycle_1.0.1    
## [46] stringi_1.7.5       yaml_2.2.1          carData_3.0-4      
## [49] MASS_7.3-54         grid_4.1.1          crayon_1.4.1       
## [52] lattice_0.20-44     haven_2.4.3         splines_4.1.1      
## [55] hms_1.1.1           pillar_1.6.4        boot_1.3-28        
## [58] estimability_1.3    ggsignif_0.6.3      reprex_2.0.1       
## [61] glue_1.4.2          evaluate_0.14       data.table_1.14.2  
## [64] modelr_0.1.8        vctrs_0.3.8         nloptr_1.2.2.2     
## [67] tzdb_0.1.2          cellranger_1.1.0    gtable_0.3.0       
## [70] assertthat_0.2.1    openxlsx_4.2.4      xfun_0.27          
## [73] xtable_1.8-4        broom_0.7.9         rstatix_0.7.0      
## [76] ellipsis_0.3.2

Themes for Plots

# Correct ggplot2 pts and fonts
ggplot_factor_line <- ggplot2::.pt*72.27/96
ggplot_factor_font = 1.0

# Create project report theme
theme_report_white  <-  theme_pubr() + theme(
  axis.line = element_line(size = 0.5 / ggplot_factor_line,color="black"),
  axis.ticks = element_line(size = 0.5 / ggplot_factor_line, color="black"),
  axis.ticks.length = unit(0.2, "cm"), 
  text = element_text(size = 8 * ggplot_factor_font,
                      family="Helvetica",
                      colour="black"),
  axis.text.x=element_text(angle = 0, 
                           vjust = 0.5,
                           size = 8 * ggplot_factor_font,
                           family = "Helvetica",
                           colour = "black"), #face="bold"
  axis.text.y = element_text(angle = 0,
                            vjust = 0.5,
                            size= 8 * ggplot_factor_font,
                            family = "Helvetica",
                            colour="black"),
  axis.title.x = element_text(angle = 0,
                            vjust = 0.5,
                            size= 8 * ggplot_factor_font,
                            family="Helvetica",
                            colour="black",
                            face="bold"),
  axis.title.y = element_text(angle = 90,
                            vjust = 0.5,
                            size= 8 * ggplot_factor_font,
                            family = "Helvetica",
                            colour = "black",
                            face = "bold") # face = "bold" 
  )

Load Data

# Demographics
df_demo <- read_excel("./Data/Data Demographics.xlsx", sheet = "Data Demo")
df_demo[c(1:5)] <- lapply(df_demo[c(1:5)], factor)

# Cognition
df_cognition <- read_excel("./Data/Data Cognition.xlsx", 
                     sheet = "Data Cognition")

# Convert var types
df_cognition[c(1:3)] <- lapply(df_cognition[c(1:3)], as.factor)

# Unify col names
names(df_cognition)[names(df_cognition) == "PVT_Accuracy"] <- "PVT_pCorr"
names(df_cognition)[names(df_cognition) == "PVT_Slowness"] <- "PVT_AvRT"
names(df_cognition)[names(df_cognition) == "BART_RiskScoreP"] <- "BART_pCorr"
names(df_cognition)[names(df_cognition) == "MP_Accuracy"] <- "MP_pCorr"
names(df_cognition)[names(df_cognition) == "LOT_Accuracy"] <- "LOT_pCorr"
names(df_cognition)[names(df_cognition) == "NBCK_Av_pCorr"] <- "NBCK_pCorr"

# Add 'z'
colnames(df_cognition)[c(4:ncol(df_cognition))] <- paste(colnames(df_cognition[c(4:ncol(df_cognition))]), "z", sep = "_")

# relabel Battery numbers
#df_cognition$time.n <- plyr::revalue(df_cognition$time.f, c("BDC-9" = "01", "BDC-7" = "02", "BDC-6" = "03", "HDT1" = "04", "HDT3" = "05", "HDT5" = "06", "HDT14" = "07", "HDT28" = "08", "HDT42" = "09", "HDT57" = "10", "R+1" = "11", "R+5" = "12", "R+12" = "13"))

#df_cognition <- df_cognition[c(1:3, 34, 4:33)]
#df_cognition$time.n <- as.character((df_cognition$time.n))
#df_cognition$time.n <- as.numeric((df_cognition$time.n))

Definition of Helper Functions

Descriptive Functions

# Function to check data completion 
func_data_completion <- function(data, ...) {
  
  if(!"time.f" %in% names(data)) {
    data$time.f <- as.factor( c( rep(1, nrow(data) ) ) )
    }
  
  if(!"cond.f" %in% names(data)) {
    data$cond.f <- as.factor( c( rep(1, nrow(data) ) ) )
    }
  
  data <- reshape2::melt( data, 
                              id.vars = names( Filter(is.factor, data) )
                              )  %>%
  
  group_by(variable) %>%
  dplyr::summarise(
      N = 24 * length( levels(time.f) ) * length( levels(cond.f) ),
      n = sum( !is.na(value) ),
      Percent = n / N * 100
      )
}

# Function to compute descriptives
func_desc <- function(data, ...) {
  df_tmp <- df_demo[ , !names(df_demo) %in% c("age","height", "weight")]
  data <- merge(df_tmp, data, by = "id") #add demo
  data <- reshape2::melt( data, 
                              id.vars = names( Filter(is.factor, data) )
                              )
  groups <- enquos(...)
  
  result <- data %>%
    
    group_by(!!!groups) %>%
    
    dplyr::summarise(
      MEAN = mean(value, na.rm = T),
      SE = sd(value, na.rm = T) / sqrt(length(value) )
      )
  return(result)
}

func_bdc_correction <- function(data, ...) {

  df_tmp <- data %>%
    subset(time.f == "BDC-6" | time.f == "BDC-7" |time.f == "BDC-9") %>%
    group_by(group_3) %>%
    dplyr::summarise(
      mean_bdc_by_group = mean(MEAN, na.rm = T)
      )
  
  df_tmp <- merge(df_tmp, data, by = "group_3")
  df_tmp %>%
    mutate(MEAN_corr = case_when(time.f == "BDC-6" | time.f == "BDC-7" | time.f == "BDC-9" ~
                              MEAN-mean_bdc_by_group,
                            TRUE ~ MEAN)) %>%
    select( c(!mean_bdc_by_group & !MEAN) ) %>%
    rename(MEAN = MEAN_corr)
}

Plot Functions

# Wrapper function for line plots, forwards to func_plot_line_facet
func_plot_line <- function(x) {ggplot(x, tmpaes_plot) + 
    theme_report_white +
    geom_line(position = pd) +
    geom_errorbar(tmpaes_errorbar,
                  width =.4, 
                  position = pd) +
    geom_point(position = pd,
               shape = 21,
               size = 2) + #shape = 21, size = 1, fill = "white"
    scale_y_continuous(expand = c(0,0)) +
    labs(y = "", x = "Time", size = 10) +
    theme(legend.title = element_blank(),
          legend.text = element_text(size = 8),
          axis.title.x = element_text(size = 10))
  }

# Function for line plots, takes 'desc' as input
func_plot_line_facet <- function(x) {func_plot_line(x) +
    facet_grid_def +
    theme(strip.placement = "outside",
          strip.background = element_rect(color = "transparent",
                                          fill = "transparent", size = 1.5,
                                          linetype = "solid"),
          strip.text = element_text(size = 10,
                                    color = "black",
                                    face = "bold"),
          panel.spacing.y = unit(2, "lines"))
}

# Function for bar plots, forwards to func_plot_bar_facet
func_plot_bar <- function(x) {ggplot(x, tmpaes_plot) + theme_report_white +
    geom_hline(yintercept = 0,
               size = 0.3,
               linetype = "solid") +
    geom_errorbar(tmpaes_errorbar,
                 width = .4,
                  position = pd,
                  color = "black") +
    geom_bar(tmpaes_geombar, # separate by category
             stat = "identity", 
             position = pd,
             width = 1.0,
             color = "black") + # make it side-by-side
    scale_y_continuous(expand = c(0,0)) +
    labs(y = "", x = "", size = 10) +
    guides(x = "none") +
    theme(legend.title = element_blank(),
          legend.text = element_text(size=8),
          axis.title.x = element_text(size = 10),
          axis.text.x = element_blank(),
          axis.ticks.x = element_blank())
}

# Function for bar plots, takes 'desc*' as input
func_plot_bar_facet <- function(x) {func_plot_bar(x) +
    facet_grid_def +
    theme(strip.placement = "outside",
          strip.background = element_rect(color = "transparent",
                                          fill = "transparent", size = 1.5,
                                          linetype = "solid"),
          strip.text = element_text(size = 10,
                                    color = "black",
                                    face = "bold"),
          panel.spacing.y = unit(2, "lines"))
}

Functions for Mixed Models

# Wrapper for appending bdc data as covariate
func_bdc_cov <- function(data, ...) {
  
  df_tmp <- subset(data, time.n > bdc )
  
  df_tmp_pre <- subset(data, time.n < (bdc+1) ) %>%
    rename_at(vars(contains("_z")), ~ str_c("pre_", .)) %>%
    select("id" | is.numeric)

  df_tmp <- merge(df_tmp_pre[-c(2)], df_tmp, by = c("id") )
  droplevels(df_tmp)
 
  df_tmp <- merge(df_demo, df_tmp, by  = c("id")) 
}

func_bdc_cov_multi_sep <- function(data, ...) {
  df_tmp <- subset(data[-c(2)], time.n > bdc )
  df_tmp_pre <- subset(data, time.n < (bdc+1) )
  df_tmp_pre <- df_tmp_pre[,!grepl("time.n", colnames(df_tmp_pre))] 
  df_tmp_pre <- 
    reshape2::melt( df_tmp_pre, 
                              id.vars = names( Filter(is.factor, df_tmp_pre) )
                              ) %>%
    group_by(id, variable) %>%
    dplyr::summarise(
      mean_bdc = mean(value, na.rm = T) 
    )
  df_tmp_pre <- reshape2::dcast( df_tmp_pre, id ~ variable, value.var = "mean_bdc") %>%
    rename_at(vars(contains("_z")), ~ str_c("pre_", .))
  df_tmp <- merge(df_tmp_pre, df_tmp, by = c("id") )
  droplevels(df_tmp)
  df_tmp <- merge(df_demo[c(1,3,5:6)], df_tmp, by  = c("id"))
  }

# Wrapper for mixed models, takes 'dvList', requires 'model_mm', e.g., model_mm <- paste0(" + sex"," + age", " + time.f",  " * group_3", " * cond.f", " + (1|id)" )
func_mm_fac_with_bdc_adj_help <- function(dvList) {
  map(dvList, function(i) {
    lmerTest::lmer(paste0(i, " ~ pre_", i, model_mm),
                   REML = T,
                   data = df_tmp)
   }) %>% setNames(dvList) #add dvList names
}

# Function for creating mixed model tables, takes 'dvList'
func_tbl_aov <- function(x) {
  tbl_mm_fac_aov <- list_df2df(tbl_mm_fac_aov)
  names(tbl_mm_fac_aov)[4] <- "df1"
  names(tbl_mm_fac_aov)[5] <- "df2"
  names(tbl_mm_fac_aov)[6] <- "F"
  names(tbl_mm_fac_aov)[7] <- "P"
  tbl_mm_fac_aov[c(5:6)] <- round(tbl_mm_fac_aov[c(5:6)], digits = 2) # Round
  tbl_mm_fac_aov[c(7)] <- round(tbl_mm_fac_aov[c(7)], digits = 3) # Round
  return(tbl_mm_fac_aov)
}

# Function emmeans, takes 'dvList'; requires to define 'model', e.g., model <- c(~ time.f * group_3 * cond.f); also requires 'func_mm_fac_with_bdc_adj_help'
func_emm_fac_with_bdc_adj <- function(x) {
  mm_fac <- func_mm_fac_with_bdc_adj_help(dvList)
  emm_fac <- lapply(mm_fac, function(mm_fac) emmeans::emmeans(mm_fac, model_emmeans, lmer.df = "satterthwaite"))
  emm_fac_ <- lapply(emm_fac, data.frame)
  table_emm_fac_df <- ldply(emm_fac_, data.frame) #create df
  table_emm_fac_df_raw <- table_emm_fac_df
  table_emm_fac_df_raw$.id <- as.factor(as.character(table_emm_fac_df_raw$.id ))
  rm(emm_fac_)
  return(table_emm_fac_df_raw)
}

# Helper function contrasts emmeans, takes 'dvList'; also requires 'func_mm_fac_with_bdc_adj_help'
func_contrasts_emm_fac_with_bdc_adj <- function(x) {
  mm_fac <- func_mm_fac_with_bdc_adj_help(dvList)
  emm_fac <- lapply(mm_fac, function(mm_fac) emmeans::emmeans(mm_fac, model_emmeans, lmer.df = "satterthwaite"))
  contrasts_mm_fac <- lapply(emm_fac, function(emm_fac) contrast(emm_fac, contrast_type,
                                                                         by = factors,
                                                                         lmer.df = "satterthwaite",
                                                                         adjust = "none"))
  contrasts_mm_fac <- lapply(contrasts_mm_fac, summary, infer=T) #add adj confidence intervals (adjusted if adjustment method is selected: https://cran.r-project.org/web/packages/emmeans/vignettes/confidence-intervals.html)
  contrasts_mm_fac_ <- lapply(contrasts_mm_fac, data.frame) #create df
  contrasts_mm_fac_df <- ldply(contrasts_mm_fac_, data.frame) #create df
  rm(contrasts_mm_fac_)
  return(contrasts_mm_fac_df)
}

# Function for creating table of contrasts emmeans, takes 'dvList', and requires 'func_contrasts_emm_fac_with_bdc_adj'
func_tbl_contr_emm_fac_with_bdc_adj <- function(dvList) {
  tmp <- func_contrasts_emm_fac_with_bdc_adj(dvList)
  tmp[ncol(tmp)] <- round(tmp[ncol(tmp)], digits = 3) #round
  tmp[ncol(tmp)-1] <- round(tmp[ncol(tmp)-1], digits = 2) #round
  tmp[ncol(tmp)-2] <- round(tmp[ncol(tmp)-2], digits = 2) #round
  tmp[ncol(tmp)-3] <- round(tmp[ncol(tmp)-3], digits = 2) #round
  tmp[ncol(tmp)-4] <- round(tmp[ncol(tmp)-4], digits = 1) #round
  tmp[ncol(tmp)-5] <- round(tmp[ncol(tmp)-5], digits = 2) #round
  tmp[ncol(tmp)-6] <- round(tmp[ncol(tmp)-6], digits = 2) #round
  tmp$CI <- paste0("(",tmp$lower.CL,"," ,tmp$upper.CL,")")
  tmp$CI <- gsub( ",", ", ", tmp$CI ) 
  tmp$CI <- paste0(tmp$estimate," " ,tmp$CI)
  tmp <- tmp[-c(ncol(tmp)-3,
                ncol(tmp)-4,
                ncol(tmp)-6,
                ncol(tmp)-7)]
  colnames(tmp)[ncol(tmp)] <- "Mean (95%CI)"
  colnames(tmp)[ncol(tmp)-1] <- "P (unadjusted)"
  colnames(tmp)[ncol(tmp)-2] <- "t"
  colnames(tmp)[2] <- "Contrast"
  return(tmp)
}

Data Completion Check

# Cognition
missing_cognition <- func_data_completion(df_cognition)
# Rename column
names(missing_cognition)[names(missing_cognition) == "variable"] <- "Variable"
names(missing_cognition)[names(missing_cognition) == "N"] <- "Expected (N)"
names(missing_cognition)[names(missing_cognition) == "n"] <- "Completed (N)"
names(missing_cognition)[names(missing_cognition) == "Percent"] <- "Completed (%)"

Table 1. Data collection summary Cognition.

## [1] 99.8

Plot Time Course adjusted for Baseline

# Plot arithmetic means
# Compute desc
desc_all <- func_desc(df_cognition, group_3, time.f, variable)

# Shift mean within groups to reflect a pre-HDT baseline performance of 0 (zero) (To reflect the analytical approach (adjusting for baseline performance).
desc_all <- func_bdc_correction(desc_all)

# Plot 
# Define aes for func_plot
tmpaes_plot <- aes(x = time.f, y = MEAN,
                   group = group_3,
                   fill = group_3,
                   color = group_3)

tmpaes_errorbar <- aes(x = time.f,
                    ymin = MEAN-SE, 
                    ymax = MEAN+SE)

# Define positioning
pd <- position_dodge(0.3)

# Speed
# Select
desc <- dplyr::filter(desc_all, grepl("AvR", variable))
desc <- droplevels(desc)
# New facet label names for variable
new_labels <- c("MP Speed", "VOLT", "NBACK", "AM", "LOT", "ERT", "MRT", "DSST", "BART", "PVT")
names(new_labels) <- c("MP_AvRT_z", "VOLT_AvRT_z", "NBCK_AvRT_z", "AM_AvRT_z", "LOT_AvRT_z", "ERT_AvRT_z", "MRT_AvRT_z", "DSST_AvRT_z", "BART_AvRT_z", "PVT_AvRT_z")

# Define labels for func_plot_line_facet
labeller_def <- labeller(variable = new_labels)

facet_grid_def <- facet_wrap(~ variable,
                             switch = "y",
                             ncol = 1,
                             labeller = labeller_def) 

# Draw plot
p.cognition_desc_AvRT <- func_plot_line_facet(desc %>%
                                                 mutate(group_3 = fct_relevel(group_3, "CTRL", "cAG", "iAG")) %>%
  mutate(time.f = fct_relevel(time.f, "BDC-9", "BDC-7", "BDC-6", "HDT1", "HDT3", "HDT5","HDT14", "HDT28", "HDT42", "HDT57", "R+1", "R+5", "R+12"))) +
  geom_hline(yintercept = 0, size = 0.3, linetype = "dotted")

# Adjust legend size
p.cognition_desc_AvRT <- p.cognition_desc_AvRT + theme(legend.key.size = unit(0.3, "cm"))

# Accuracy
# Select
desc <- dplyr::filter(desc_all, grepl("pCorr",variable))
desc <- droplevels(desc)

# New facet label names for variable
new_labels <- c("MP", "VOLT", "NBACK", "AM", "LOT", 
                      "ERT", "MRT", "DSST", "BART", "PVT")
names(new_labels) <- c("MP_pCorr_z", "VOLT_pCorr_z", "NBCK_Av_pCorr_z", "AM_pCorr_z", "LOT_pCorr_z", "ERT_pCorr_z", "MRT_pCorr_z", "DSST_pCorr_z", "BART_pCorr_z", "PVT_pCorr_z")

# Define labels for func_plot_line_facet
labeller_def <- labeller(variable = new_labels)

facet_grid_def <- facet_wrap(~ variable,
                             switch = "y",
                             ncol = 1,
                             labeller = labeller_def) 

# Draw plot
p.cognition_desc_pCorr <- func_plot_line_facet(desc %>%
                                                 mutate(group_3 = fct_relevel(group_3, "CTRL", "cAG", "iAG")) %>%
  mutate(time.f = fct_relevel(time.f, "BDC-9", "BDC-7", "BDC-6", "HDT1", "HDT3", "HDT5","HDT14", "HDT28", "HDT42", "HDT57", "R+1", "R+5", "R+12"))) +
  geom_hline(yintercept = 0, size = 0.3, linetype = "dotted")

# Adjust legend size
p.cognition_desc_pCorr <- p.cognition_desc_pCorr + theme(legend.key.size = unit(0.3, "cm"))

# Efficiency
# Select
desc <- dplyr::filter(desc_all, grepl("Eff",variable))
desc <- subset(desc, variable != "BART_Eff_z")
desc <- droplevels(desc)

# New facet label names for variable
new_labels <- c("MP", "VOLT", "NBACK", "AM", "LOT", 
                      "ERT", "MRT", "DSST" , "PVT")
names(new_labels) <- c("MP_Eff_z", "VOLT_Eff_z", "NBCK_Eff_z", "AM_Eff_z", "LOT_Eff_z", "ERT_Eff_z", "MRT_Eff_z", "DSST_Eff_z", "PVT_Eff_z")

# Define labels for func_plot_line_facet
labeller_def <- labeller(variable = new_labels)

facet_grid_def <- facet_wrap(~ variable,
                             switch = "y",
                             ncol = 1,
                             labeller = labeller_def) 

# Draw plot
p.cognition_desc_eff <- func_plot_line_facet(desc %>%
                                                 mutate(group_3 = fct_relevel(group_3, "CTRL", "cAG", "iAG")) %>%
  mutate(time.f = fct_relevel(time.f, "BDC-9", "BDC-7", "BDC-6", "HDT1", "HDT3", "HDT5","HDT14", "HDT28", "HDT42", "HDT57", "R+1", "R+5", "R+12"))) +
  geom_hline(yintercept = 0, size = 0.3, linetype = "dotted")

# Adjust legend size
p.cognition_desc_eff <- p.cognition_desc_eff + theme(legend.key.size = unit(0.3, "cm"))

Mixed Model Time x Group Interaction

# Create baseline data
# relabel time.f
df_cognition$time.n <- plyr::revalue(df_cognition$time.f, c("BDC-9" = "01", "BDC-7" = "02", "BDC-6" = "03", "HDT1" = "04", "HDT3" = "05", "HDT5" = "06", "HDT14" = "07", "HDT28" = "08", "HDT42" = "09", "HDT57" = "10", "R+1" = "11", "R+5" = "12", "R+12" = "13"))
df_cognition <- df_cognition[c(1:3, 34, 4:33)]
df_cognition$time.n <- as.character((df_cognition$time.n))
df_cognition$time.n <- as.numeric((df_cognition$time.n))

# Define BDC sessions (number indicates last bdc sessions based on time.n)
bdc <- 3

# Create df with bdc as covariate
df_tmp <- func_bdc_cov_multi_sep(df_cognition, id)

# Select dep vars
dvList <- colnames(df_tmp)[c(37:66)]

# Define model
model_mm <- paste0(" + sex"," + age", " + time.f",  " * group_3", " + (1|id)" )

# Run mixed models
tbl_mm_fac_aov <- lapply(func_mm_fac_with_bdc_adj_help(dvList), anova, simply = FALSE, USE.NAMES = TRUE) 

# ANOVA table 
tbl_mm_fac_aov_df <- func_tbl_aov(dvList)

# Sort by variable
tbl_mm_fac_aov_df_AvRT <- dplyr::filter(tbl_mm_fac_aov_df, grepl("AvR", X1))
tbl_mm_fac_aov_df_pCorr <- dplyr::filter(tbl_mm_fac_aov_df, grepl("pCorr", X1))
tbl_mm_fac_aov_df_eff <- dplyr::filter(tbl_mm_fac_aov_df, grepl("Eff", X1))

# Helper functions for parameter definition in mixed model tables, takes 'dvList'
Parameter_func <- function(dvList){
  rep(c(
  "Baseline",
  "Sex",
  "Age",
  "Time",
  "Group",
  "Time x Group"),
                 length(dvList)/3)
}

# Create parameters
Parameter <- Parameter_func(dvList)
blanks <- rep("", length(Parameter) / ( length(dvList) / 3) - 1)
Variable <-  c("MP", blanks,"VOLT", blanks, "NBCK", blanks, "AM", blanks, "LOT", blanks, "ERT", blanks, "MRT", blanks, "DSST", blanks, "BART", blanks, "PVT", blanks)

# Add to parameters to tables
tbl_mm_fac_aov_df_AvRT <- cbind(Variable, Parameter, tbl_mm_fac_aov_df_AvRT[c(4:7)])
tbl_mm_fac_aov_df_pCorr <- cbind(Variable, Parameter, tbl_mm_fac_aov_df_pCorr[c(4:7)])
tbl_mm_fac_aov_df_eff <- cbind(Variable, Parameter, tbl_mm_fac_aov_df_eff[c(4:7)])

# Drop BART from Eff
tbl_mm_fac_aov_df_eff <- tbl_mm_fac_aov_df_eff[-c( 49:54),]

Table 2. Mixed Model Summary for Speed.

Table 3. Mixed Model Summary for Accuracy.

Table 4. Mixed Model Summary for Efficiency.

Contrasts for Comparing Groups at Each Time Point During Bed Rest

# Contrasts
# Define model for contrasts
model_emmeans <- ~ time.f * group_3

# Define contrasts
contrast_type <- "pairwise"
factors <- c("time.f")

# Run contrasts
tbl_contr_emm_fac_df <- func_tbl_contr_emm_fac_with_bdc_adj(dvList %>% 
                                       mutate(time.f = fct_relevel(time.f, "HDT1", "HDT3", "HDT5", "HDT14", "HDT28", "HDT42","HDT57", "R+1", "R+5", "R+12"))) 

# Add variables
# AvRT
tbl_contr_emm_fac_df_AvRT <- dplyr::filter(tbl_contr_emm_fac_df, grepl("AvR", .id))
#tbl_contr_emm_fac_df_AvRT <- tbl_contr_emm_fac_df_AvRT[c(8:9, 2, 7, 4:6)]

# pCorr
tbl_contr_emm_fac_df_pCorr <- dplyr::filter(tbl_contr_emm_fac_df, grepl("pCorr", .id))
#tbl_contr_emm_fac_df_pCorr <- tbl_contr_emm_fac_df_pCorr[c(8:9, 2, 7, 4:6)]

# Eff
tbl_contr_emm_fac_df_eff <- dplyr::filter(tbl_contr_emm_fac_df, grepl("Eff", .id))
#tbl_contr_emm_fac_df_eff <- tbl_contr_emm_fac_df_eff[c(8:9, 2, 7, 4:6)]

# Create vars to make nice table
tmp_n_dvs <- length(levels(tbl_contr_emm_fac_df$Contrast)) * 
  length(levels(tbl_contr_emm_fac_df$time.f)) -1

tmp_n_time <- length(levels(tbl_contr_emm_fac_df$Contrast)) - 1

Variable <- c("MP", rep("", tmp_n_dvs),
                                   "VOLT", rep("", tmp_n_dvs),
                                   "NBCK", rep("", tmp_n_dvs),
                                   "AM", rep("", tmp_n_dvs),
                                   "LOT", rep("", tmp_n_dvs),
                                   "ERT", rep("", tmp_n_dvs),
                                   "MRT", rep("", tmp_n_dvs),
                                   "DSST", rep("", tmp_n_dvs),
                                   "BART", rep("", tmp_n_dvs),
                                   "PVT", rep("", tmp_n_dvs))

Time <- rep ( c("HDT1", rep("", tmp_n_time),
                                     "HDT3", rep("", tmp_n_time),
                                     "HDT5", rep("", tmp_n_time),
                                     "HDT14", rep("", tmp_n_time),
                                     "HDT28", rep("", tmp_n_time),
                                     "HDT42", rep("", tmp_n_time),
                                     "HDT57", rep("", tmp_n_time),
                                     "R+1", rep("", tmp_n_time),
                                     "R+5", rep("", tmp_n_time),
                                     "R+12", rep("", tmp_n_time)),
                                   length(dvList)/3 )

# Add to parameters to tables
tbl_contr_emm_fac_df_AvRT <- cbind(Variable, Time, tbl_contr_emm_fac_df_AvRT[c(2, 7, 4:6)])

tbl_contr_emm_fac_df_pCorr <- cbind(Variable, Time, tbl_contr_emm_fac_df_pCorr[c(2, 7, 4:6)])

tbl_contr_emm_fac_df_eff <- cbind(Variable, Time, tbl_contr_emm_fac_df_eff)
tbl_contr_emm_fac_df_eff <- dplyr::filter(tbl_contr_emm_fac_df_eff, !grepl("BART_Eff_z", .id)) # Drop BART from Eff
tbl_contr_emm_fac_df_eff <- tbl_contr_emm_fac_df_eff[c(1,2,4,5,9,6:8)]

Table 5. Mixed Model Group Contrasts for Speed.

Table 6. Mixed Model Group Contrasts for Accuracy.

Table 7. Mixed Model Group Contrasts for Efficiency

Plot Estimated Marignal Means To Show Time Course for Each Group

# Define model for emmeans
model_emmeans <- ~ time.f * group_3

# Run emmeans
emm_fac <- func_emm_fac_with_bdc_adj(dvList) 

# Plot emmeans 
# Define aes for func_plot
tmpaes_plot <- aes(x = time.f, y = emmean,
                   group = group_3,
                   fill = group_3,
                   color = group_3)

tmpaes_errorbar <- aes(x = time.f,
                    ymin = emmean-SE, 
                    ymax = emmean+SE)

# Draw plot
pd <- position_dodge(0.3)

# Accuracy
# Select
emm_fac_pCorr <- dplyr::filter(emm_fac, grepl("pCorr", .id))
emm_fac_pCorr <- droplevels(emm_fac_pCorr)

# New facet label names for variable
new_labels <- c("MP", "VOLT", "NBACK", "AM", "LOT", 
                      "ERT", "MRT", "DSST", "BART", "PVT")
names(new_labels) <- c("MP_pCorr_z", "VOLT_pCorr_z", "NBCK_Av_pCorr_z", "AM_pCorr_z", "LOT_pCorr_z", "ERT_pCorr_z", "MRT_pCorr_z", "DSST_pCorr_z", "BART_pCorr_z", "PVT_pCorr_z")

# Define labels for func_plot_line_facet
labeller_def <- labeller(.id = new_labels)

# Define grid
facet_grid_def <- facet_wrap(~ .id,
                             switch = "y",
                             ncol = 1,
                             labeller = labeller_def) 

# Plot
p.cognition_emmeans_pCorr <- func_plot_line_facet(emm_fac_pCorr %>%
                                                  mutate(group_3 = fct_relevel(group_3, "CTRL", "cAG", "iAG"),
                                                         .id = fct_relevel(.id, "MP_pCorr_z", "VOLT_pCorr_z", "NBCK_Av_pCorr_z","AM_pCorr_z","LOT_pCorr_z","ERT_pCorr_z","MRT_pCorr_z", "DSST_pCorr_z", "BART_pCorr_z", "PVT_pCorr_z")))+ geom_hline(yintercept = 0, size = 0.3, linetype = "dotted") + theme(legend.key.size = unit(0.3, "cm"))

# Speed
# Select
emm_fac_AvRT <- dplyr::filter(emm_fac, grepl("AvR", .id))
emm_fac_AvRT <- droplevels(emm_fac_AvRT)

# New facet label names for variable
new_labels <- c("MP Speed", "VOLT", "NBACK", "AM", "LOT", 
                      "ERT", "MRT", "DSST", "BART", "PVT")
names(new_labels) <- c("MP_AvRT_z", "VOLT_AvRT_z", "NBCK_AvRT_z", "AM_AvRT_z", "LOT_AvRT_z", "ERT_AvRT_z", "MRT_AvRT_z", "DSST_AvRT_z", "BART_AvRT_z", "PVT_AvRT_z")

# Define labels for func_plot_line_facet
labeller_def <- labeller(.id = new_labels)

# Define grid
facet_grid_def <- facet_wrap(~ .id,
                             switch = "y",
                             ncol = 1,
                             labeller = labeller_def) 

# Plot
p.cognition_emmeans_AvRT <- func_plot_line_facet(emm_fac_AvRT %>%
                                                  mutate(group_3 = fct_relevel(group_3, "CTRL", "cAG", "iAG"),
                                                         .id = fct_relevel(.id, "MP_AvRT_z", "VOLT_AvRT_z", "NBCK_AvRT_z","AM_AvRT_z","LOT_AvRT_z","ERT_AvRT_z","MRT_AvRT_z", "DSST_AvRT_z", "BART_AvRT_z", "PVT_AvRT_z"))) + geom_hline(yintercept = 0, size = 0.3, linetype = "dotted") + theme(legend.key.size = unit(0.3, "cm"))

# Efficiency
# Select
emm_fac_eff <- dplyr::filter(emm_fac, grepl("Eff", .id))
emm_fac_eff <- droplevels(emm_fac_eff)
# Drop BART and PVT from Eff
emm_fac_eff <- emm_fac_eff %>% filter(!grepl("BART_Eff_z", .id))

# New facet label names for variable
new_labels <- c("MP", "VOLT", "NBACK", "AM", "LOT", 
                      "ERT", "MRT", "DSST", "PVT")
names(new_labels) <- c("MP_Eff_z", "VOLT_Eff_z", "NBCK_Eff_z", "AM_Eff_z", "LOT_Eff_z", "ERT_Eff_z", "MRT_Eff_z", "DSST_Eff_z", "PVT_Eff_z")

# Define labels for func_plot_line_facet
labeller_def <- labeller(.id = new_labels)

# Define grid
facet_grid_def <- facet_wrap(~ .id,
                             switch = "y",
                             ncol = 1,
                             labeller = labeller_def) 

# Plot
p.cognition_emmeans_eff <- func_plot_line_facet(emm_fac_eff %>%
                                                  mutate(group_3 = fct_relevel(group_3, "CTRL", "cAG", "iAG"),
                                                         .id = fct_relevel(.id, "MP_Eff_z", "VOLT_Eff_z", "NBCK_Eff_z","AM_Eff_z","LOT_Eff_z","ERT_Eff_z","MRT_Eff_z", "DSST_Eff_z", "PVT_Eff_z"))) + geom_hline(yintercept = 0, size = 0.3, linetype = "dotted") + theme(legend.key.size = unit(0.3, "cm"))

Mixed Model for Main Effect of Bed Rest (excluding BDC)

# Select HDBR data
# Drop Recovery
df_tmp <- subset(df_tmp, time.n < 11)
df_tmp <- droplevels(df_tmp)

# Select dep vars
dvList <- colnames(df_tmp)[c(37:66)]

# Define model
model_mm <- paste0(" + sex"," + age", " + group_3", " + (1|id)" )

# Run mixed models
tbl_mm_fac_aov <- lapply(func_mm_fac_with_bdc_adj_help(dvList), anova, simply = FALSE, USE.NAMES = TRUE) 

# ANOVA table
tbl_mm_fac_aov_df <- func_tbl_aov(dvList)

# Helper functions for parameter definition in mixed model tables, takes 'dvList'
Parameter_func <- function(dvList){
  rep(c(
  "Baseline",
  "Sex",
  "Age",
  "Group"),
                 length(dvList)/3)
}

# Add variables
tbl_mm_fac_aov_df_AvRT <- dplyr::filter(tbl_mm_fac_aov_df, grepl("AvR", X1))

tbl_mm_fac_aov_df_pCorr <- dplyr::filter(tbl_mm_fac_aov_df, grepl("pCorr", X1))

tbl_mm_fac_aov_df_eff <- dplyr::filter(tbl_mm_fac_aov_df, grepl("Eff", X1))


# Create parameters
Parameter <- Parameter_func(dvList)
blanks <- rep("", length(Parameter) / ( length(dvList)/3 ) - 1)
Variable <- c("MP", blanks,"VOLT", blanks, "NBCK", blanks, "AM", blanks, "LOT", blanks, "ERT", blanks, "MRT", blanks, "DSST", blanks, "BART", blanks, "PVT", blanks)

# Add to parameters to tables
tbl_mm_fac_aov_df_AvRT <- cbind(Variable, Parameter, tbl_mm_fac_aov_df_AvRT[c(4:7)])
tbl_mm_fac_aov_df_pCorr <- cbind(Variable, Parameter, tbl_mm_fac_aov_df_pCorr[c(4:7)])

tbl_mm_fac_aov_df_eff <- cbind(Variable, Parameter, tbl_mm_fac_aov_df_eff)
tbl_mm_fac_aov_df_eff <- dplyr::filter(tbl_mm_fac_aov_df_eff, !grepl("BART_Eff_z", X1)) # Drop BART from Eff
tbl_mm_fac_aov_df_eff <- tbl_mm_fac_aov_df_eff[-c(3:5)]

Table 8. Mixed Model Summary testing the Main Effect of Group for Speed.

Table 9. Mixed Model Summary testing the Main Effect of Group for Accuracy.

Table 10. Mixed Model Summary testing the Main Effect of Group for Efficiency.

Contrasts for Comparing Groups

# Define model for emmeans
model_emmeans <- ~ group_3

# Run emmeans
emm_fac <- func_emm_fac_with_bdc_adj(dvList)

# Contrasts
# Define model for contrasts
model_emmeans <- ~ group_3

# Define contrasts
contrast_type <- "pairwise"
factors <- c()

# Run contrasts
tbl_contr_emm_fac_df <- func_tbl_contr_emm_fac_with_bdc_adj(dvList)

# Add variables
# AvRT
tbl_contr_emm_fac_df_AvRT <- dplyr::filter(tbl_contr_emm_fac_df, grepl("AvR", .id))
#tbl_contr_emm_fac_df_AvRT <- tbl_contr_emm_fac_df_AvRT[c(8:9, 2, 7, 4:6)]

# pCorr
tbl_contr_emm_fac_df_pCorr <- dplyr::filter(tbl_contr_emm_fac_df, grepl("pCorr", .id))
#tbl_contr_emm_fac_df_pCorr <- tbl_contr_emm_fac_df_pCorr[c(8:9, 2, 7, 4:6)]

# Eff
tbl_contr_emm_fac_df_eff <- dplyr::filter(tbl_contr_emm_fac_df, grepl("Eff", .id))
#tbl_contr_emm_fac_df_eff <- tbl_contr_emm_fac_df_eff[c(8:9, 2, 7, 4:6)]


# Create vars to make nice table
Variable <- c("MP", rep("", 2),
                                   "VOLT", rep("", 2),
                                   "NBCK", rep("", 2),
                                   "AM", rep("", 2),
                                   "LOT", rep("", 2),
                                   "ERT", rep("", 2),
                                   "MRT", rep("", 2),
                                   "DSST", rep("", 2),
                                   "BART", rep("", 2),
                                   "PVT", rep("", 2))

# Add to parameters to tables
tbl_contr_emm_fac_df_AvRT <- cbind(Variable, tbl_contr_emm_fac_df_AvRT[c(2, 6, 3:5)])

tbl_contr_emm_fac_df_pCorr <- cbind(Variable, tbl_contr_emm_fac_df_pCorr[c(2, 6, 3:5)])

tbl_contr_emm_fac_df_eff <- cbind(Variable, tbl_contr_emm_fac_df_eff[c(2, 6, 3:5)])

# Drop BART from Eff
tbl_contr_emm_fac_df_eff <- tbl_contr_emm_fac_df_eff[-c( 25:27),]

Table 11. Mixed Model Contrasts Comparing Groups for Speed.

Table 12. Mixed Model Contrasts Comparing Groups for Accuracy.

Table 13. Mixed Model Contrasts Comparing Groups for Efficiency.

Plot Estimated Marignal Means for Comparing Groups

# Define aes for func_plot
tmpaes_plot <- aes(x = contrast,
                   y = emmean,
                   group = group_3,
                   fill = group_3,
                   color = group_3)

tmpaes_errorbar <- aes(x = group_3,
                    ymin = emmean-SE, 
                    ymax = emmean+SE)

tmpaes_geombar <- aes(x = group_3,
                    y = emmean, 
                    fill=group_3)

# Positioning
pd <- position_dodge(0.99)

# AvRT
emm_fac_AvRT <- dplyr::filter(emm_fac, grepl("AvR", .id))
emm_fac_AvRT <- droplevels(emm_fac_AvRT)

# New facet label names for variable
new_labels <- c("MP", "VOLT", "NBACK", "AM", "LOT", 
                      "ERT", "MRT", "DSST", "BART", "PVT")
names(new_labels) <- c("MP_AvRT_z", "VOLT_AvRT_z", "NBCK_AvRT_z", "AM_AvRT_z", "LOT_AvRT_z", "ERT_AvRT_z", "MRT_AvRT_z", "DSST_AvRT_z", "BART_AvRT_z", "PVT_AvRT_z")

# Define labels for func_plot_line_facet
labeller_def <- labeller(.id = new_labels)

# Define grid
facet_grid_def <- facet_wrap(~ .id,
                             switch = "y",
                             ncol = 1,
                             labeller = labeller_def) 

# Draw Plot 
p.cognition_emmeans_AvRT_main <- func_plot_bar_facet(emm_fac_AvRT %>%
                                                  mutate(group_3 = fct_relevel(group_3, "CTRL", "cAG", "iAG"),
                                                         .id = fct_relevel(.id, "MP_AvRT_z", "VOLT_AvRT_z", "NBCK_AvRT_z","AM_AvRT_z","LOT_AvRT_z","ERT_AvRT_z","MRT_AvRT_z", "DSST_AvRT_z", "BART_AvRT_z", "PVT_AvRT_z"))) + theme(legend.key.size = unit(0.3, "cm"))


# Accuracy
# Select
emm_fac_pCorr <- dplyr::filter(emm_fac, grepl("pCorr", .id))
emm_fac_pCorr <- droplevels(emm_fac_pCorr)

# New facet label names for variable
new_labels <- c("MP", "VOLT", "NBACK", "AM", "LOT", 
                      "ERT", "MRT", "DSST", "BART", "PVT")
names(new_labels) <- c("MP_pCorr_z", "VOLT_pCorr_z", "NBCK_Av_pCorr_z", "AM_pCorr_z", "LOT_pCorr_z", "ERT_pCorr_z", "MRT_pCorr_z", "DSST_pCorr_z", "BART_pCorr_z", "PVT_pCorr_z")

# Define labels for func_plot_line_facet
labeller_def <- labeller(.id = new_labels)

# Define grid
facet_grid_def <- facet_wrap(~ .id,
                             switch = "y",
                             ncol = 1,
                             labeller = labeller_def) 

# Draw Plot 
p.cognition_emmeans_pCorr_main <- func_plot_bar_facet(emm_fac_pCorr %>%
                                                  mutate(group_3 = fct_relevel(group_3, "CTRL", "cAG", "iAG"),
                                                         .id = fct_relevel(.id, "MP_pCorr_z", "VOLT_pCorr_z", "NBCK_Av_pCorr_z","AM_pCorr_z","LOT_pCorr_z","ERT_pCorr_z","MRT_pCorr_z", "DSST_pCorr_z", "BART_pCorr_z", "PVT_pCorr_z"))) + theme(legend.key.size = unit(0.3, "cm"))


# Efficiency
# Select
emm_fac_eff <- dplyr::filter(emm_fac, grepl("Eff", .id))
emm_fac_eff <- droplevels(emm_fac_eff)
# Drop BART from Eff
emm_fac_eff <- emm_fac_eff[-c( 25:27),]

# New facet label names for variable
new_labels <- c("MP", "VOLT", "NBACK", "AM", "LOT", 
                      "ERT", "MRT", "DSST", "PVT")
names(new_labels) <- c("MP_Eff_z", "VOLT_Eff_z", "NBCK_Eff_z", "AM_Eff_z", "LOT_Eff_z", "ERT_Eff_z", "MRT_Eff_z", "DSST_Eff_z", "PVT_Eff_z")

# Define labels for func_plot_line_facet
labeller_def <- labeller(.id = new_labels)

# Define grid
facet_grid_def <- facet_wrap(~ .id,
                             switch = "y",
                             ncol = 1,
                             labeller = labeller_def) 

# Draw Plot 
p.cognition_emmeans_eff_main <- func_plot_bar_facet(emm_fac_eff %>%
                                                  mutate(group_3 = fct_relevel(group_3, "CTRL", "cAG", "iAG"),
                                                         .id = fct_relevel(.id, "MP_Eff_z", "VOLT_Eff_z", "NBCK_Eff_z","AM_Eff_z","LOT_Eff_z","ERT_Eff_z","MRT_Eff_z", "DSST_Eff_z", "PVT_Eff_z"))) + theme(legend.key.size = unit(0.3, "cm"))

Create Combined plot Speed

# Combine plot 
p.cognition_comb_speed <- plot_grid(p.cognition_desc_AvRT,
                                  p.cognition_emmeans_AvRT_main + rremove("x.text") + theme( strip.text.y = element_blank() ),
                                  labels = c("a", "b"),
                                  ncol = 2, nrow = 1, rel_widths = c(2.2, 0.8)
                                  )
# Adjust legend size
p.cognition_comb_speed <- p.cognition_comb_speed + theme(legend.key.size = unit(0.3, "cm"))

# Save plot as pdf
pdf("./Figures/Cognition/p.cognition_comb_speed.pdf",
    width = 32/2.54, 
    height = 22/2.54,
    useDingbats=FALSE) #in cm
print(p.cognition_comb_speed)
dev.off()

## quartz_off_screen 
##                 2

Figure 1. (a) Time course of speed before, during, and after bed rest. Means within groups are shifted to reflect a pre-HDT baseline performance of 0 (zero). (b) Main effect of bed rest. Data means and SE are adjusted for baseline.

## Create Combined plot Accuracy

# Combine plot 
p.cognition_comb_acc <- plot_grid(p.cognition_desc_pCorr,
                                  p.cognition_emmeans_pCorr_main + rremove("x.text") + theme( strip.text.y = element_blank() ),
                                  labels = c("a", "b"),
                                  ncol = 2, nrow = 1, rel_widths = c(2.2, 0.8)
                                  )

# Adjust legend size
p.cognition_comb_acc <- p.cognition_comb_acc + theme( legend.key.size = unit(0.3, "cm") )

# Save plot as pdf
pdf("./Figures/Cognition/p.cognition_comb_acc.pdf",
    width = 22/2.54, 
    height = 32/2.54,
    useDingbats=FALSE) #in cm
print(p.cognition_comb_acc)
dev.off()

## quartz_off_screen 
##                 2

Figure 2. (a) Time course of accuracy before, during, and after bed rest. Means within groups are shifted to reflect a pre-HDT baseline performance of 0 (zero). (b) Main effect of bed rest. Data means and SE are adjusted for baseline.

Create Combined plot Efficiency

# Combine plot 
p.cognition_eff_comb <- plot_grid(p.cognition_desc_eff, p.cognition_emmeans_eff_main + rremove("x.text") + theme( strip.text.y = element_blank() ),
                                  labels = c("a", "b"),
                                  ncol = 2, nrow = 1, rel_widths = c(2.2, 0.8)
                                  )

# Adjust legend size
p.cognition_eff_comb <- p.cognition_eff_comb + theme(legend.key.size = unit(0.3, "cm"))

# Save plot as pdf
pdf("./Figures/Cognition/p.cognition_eff_comb.pdf",
    width = 22/2.54, 
    height = 30/2.54,
    useDingbats=FALSE) #in cm
print(p.cognition_eff_comb)
dev.off()

## quartz_off_screen 
##                 2

Figure 3. (a) Time course of efficiency before, during, and after bed rest. Means within groups are shifted to reflect a pre-HDT baseline performance of 0 (zero). (b) Main effect of bed rest. Data means and SE are adjusted for baseline.